function X = tensorfun(op,A,B)
%TENSORFUN
%
%   X = tensorfun(operation,A,B)
%
%Generalization of KRON (as modified, e.g., by Laurent Sorber, Bruno Luong, and others)
%
%Creates matrix X consisting of blocks X_ij=bsxfun(operation,a(i,j),B).
%Note that "operation" must preserve the size of B.
[I J] = size(A);
[K L] = size(B);
if ~issparse(A) && ~issparse(B)
    A = reshape(A,[1 I 1 J]);
    B = reshape(B,[K 1 L 1]);
    X = reshape(bsxfun(op,A,B),[I*K J*L]);
else
      [ia,ja,sa] = find(A); ia=ia(:); ja=ja(:); sa=sa(:);
      [ib,jb,sb] = find(B); ib=ib(:); jb=jb(:); sb=sb(:);
      ix = bsxfun(op,K*(ia-1).',ib);
      jx = bsxfun(op,L*(ja-1).',jb);
      X = sparse(ix,jx,bsxfun(op,sb,sa.'),I*K,J*L);
end